import cv2
import numpy as np
import os
from scipy.optimize import curve_fit
from scipy.signal import savgol_filter
from scipy.fft import fft, fftfreq
from scipy.interpolate import splprep, splev, interp1d
from skimage.measure import label, regionprops
from skimage.morphology import skeletonize
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl
import math
import warnings
import traceback
from typing import List, Tuple, Dict, Any, Optional

warnings.filterwarnings('ignore')

# ==================== 全局配置 ====================
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans', 'Arial Unicode MS', 'Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False

# 根据文档1定义的钻孔参数
BOREHOLE_DIAMETER = 30.0  # mm
BOREHOLE_CIRCUMFERENCE = np.pi * BOREHOLE_DIAMETER  # ≈94.25 mm
BOREHOLE_DEPTH = 500.0  # mm (单孔深度)

# 根据文档1定义的DPI参数 - 仅使用单孔DPI
SINGLE_BOREHOLE_DPI = 65.73  # 单孔图像DPI
MM_PER_INCH = 25.4  # 1英寸=25.4毫米

# 输入输出路径
INPUT_DIR = "dataset/result/"
OUTPUT_DIR = "results/"
os.makedirs(OUTPUT_DIR, exist_ok=True)

# 采样方法配置
SAMPLING_METHODS = ['equidistant', 'arc_length', 'random', 'curvature_based']  # 支持的采样方法列表
DEFAULT_NUM_POINTS = 100  # 默认采样点数量，控制变量


# ==================== 核心函数定义 ====================

def sine_func(x: np.ndarray, R: float, P: float, beta: float, C: float) -> np.ndarray:
    """正弦函数模型：y = R * sin(2πx/P + beta) + C"""
    return R * np.sin(2 * np.pi * x / P + beta) + C


def pixel_to_mm_coordinates(x_pixel: np.ndarray, y_pixel: np.ndarray,
                            image_width: int, image_height: int) -> Tuple[np.ndarray, np.ndarray]:
    """
    将像素坐标转换为毫米坐标（使用固定单孔DPI）

    参数:
        x_pixel: x方向像素坐标
        y_pixel: y方向像素坐标
        image_width: 图像宽度（像素）
        image_height: 图像高度（像素）

    返回:
        Tuple[np.ndarray, np.ndarray]: (x_mm, y_mm) 毫米坐标
    """
    # 计算毫米与像素的转换系数
    mm_per_pixel_x = MM_PER_INCH / SINGLE_BOREHOLE_DPI
    mm_per_pixel_y = MM_PER_INCH / SINGLE_BOREHOLE_DPI

    # 转换为毫米坐标
    x_mm = x_pixel * mm_per_pixel_x
    y_mm = y_pixel * mm_per_pixel_y

    return x_mm, y_mm


def enhance_image(image: np.ndarray) -> np.ndarray:
    """增强图像质量，减少噪声并连接断裂的裂隙"""
    try:
        filtered = cv2.medianBlur(image, 3)
        kernel = np.ones((3, 3), np.uint8)
        return cv2.morphologyEx(filtered, cv2.MORPH_CLOSE, kernel)
    except Exception:
        return image


def validate_fracture(coords: np.ndarray, min_length_mm=30, min_x_span_mm=30.0,
                      min_y_std_mm=2.0) -> Tuple[bool, str]:
    """验证识别出的区域是否为有效的正弦状裂隙（使用毫米单位）"""
    if len(coords) == 0:
        return False, "坐标数组为空"

    x_data, y_data = coords[:, 0], coords[:, 1]

    # 检查物理尺寸（毫米）
    if len(x_data) < 10:  # 最小像素数量
        return False, f"像素数量不足({len(x_data)})"

    # 检查x方向跨度（毫米）
    x_span = np.max(x_data) - np.min(x_data)
    if x_span < min_x_span_mm:
        return False, f"x方向跨度不足({x_span:.2f}mm < {min_x_span_mm}mm)"

    # 检查y方向变化（毫米）
    y_std = np.std(y_data)
    if y_std < min_y_std_mm:
        return False, f"y方向变化太小({y_std:.2f}mm < {min_y_std_mm}mm)"

    return True, "验证通过"


def calculate_r_squared(y_true: np.ndarray, y_pred: np.ndarray) -> float:
    """计算决定系数R²"""
    ss_res = np.sum((y_true - y_pred) ** 2)
    ss_tot = np.sum((y_true - np.mean(y_true)) ** 2)
    return 1 - (ss_res / ss_tot) if ss_tot != 0 else 0.0


def robust_curve_fit(x_data: np.ndarray, y_data: np.ndarray,
                     R_init: float, P_init: float, beta_init: float, C_init: float) -> tuple:
    """鲁棒的曲线拟合 - 多策略优化以避免局部最优解"""
    best_params = None
    best_r2 = -np.inf

    # 更合理的参数边界
    R_min = max(0.1, R_init * 0.1)
    R_max = min(100.0, R_init * 5.0)

    P_min = max(50.0, P_init * 0.5)
    P_max = min(150.0, P_init * 2.0)

    C_min = np.min(y_data) - 20
    C_max = np.max(y_data) + 20

    bounds = (
        [R_min, P_min, -2 * np.pi, C_min],
        [R_max, P_max, 2 * np.pi, C_max]
    )

    # 确保初始值在边界内
    R_init_clipped = np.clip(R_init, R_min, R_max)
    P_init_clipped = np.clip(P_init, P_min, P_max)
    C_init_clipped = np.clip(C_init, C_min, C_max)

    p0 = [R_init_clipped, P_init_clipped, beta_init, C_init_clipped]

    try:
        popt, _ = curve_fit(sine_func, x_data, y_data, p0=p0, bounds=bounds, maxfev=10000)

        R, P, beta, C = popt
        beta = (beta + np.pi) % (2 * np.pi) - np.pi

        y_pred = sine_func(x_data, R, P, beta, C)
        r_squared = calculate_r_squared(y_data, y_pred)

        return (R, P, beta, C, r_squared)
    except Exception as e:
        print(f"曲线拟合失败: {str(e)}")
        return (R_init, P_init, beta_init, C_init, 0.0)


def extract_contour_skeleton(region_mask: np.ndarray) -> np.ndarray:
    """
    使用骨架化算法提取裂隙区域的轮廓线。
    参数:
        region_mask: 二值掩码图像，其中裂隙区域为1，背景为0。
    返回:
        contour_points: 轮廓线离散坐标点数组，形状为(N, 2)，其中每行是(x, y)像素坐标。
    """
    try:
        # 确保输入是二值图像
        if region_mask.dtype != np.uint8:
            region_mask = region_mask.astype(np.uint8)

        # 应用骨架化算法
        skeleton = skeletonize(region_mask > 0)

        # 获取骨架点的坐标
        y_coords, x_coords = np.where(skeleton)
        contour_points = np.column_stack((x_coords, y_coords))

        # 如果点太少，返回空数组
        if len(contour_points) < 2:
            return np.array([])

        return contour_points
    except Exception as e:
        print(f"骨架化提取错误: {str(e)}")
        return np.array([])


def calculate_jrc(x_points: np.ndarray, y_points: np.ndarray,
                  sampling_method: str = 'equidistant', num_points: int = DEFAULT_NUM_POINTS) -> Tuple[float, float]:
    """
    计算JRC值基于离散坐标点，支持四种采样方法。
    参数:
        x_points: x坐标数组。
        y_points: y坐标数组。
        sampling_method: 采样方法，可选 'equidistant' (等间距), 'arc_length' (等弧长), 'random' (随机), 'curvature_based' (基于曲率)。
        num_points: 采样点数量。
    返回:
        jrc: JRC值。
        z2: Z2值。
    """
    try:
        # 预处理坐标点：排序并去重
        sorted_indices = np.argsort(x_points)
        x_sorted = x_points[sorted_indices]
        y_sorted = y_points[sorted_indices]

        # 去除重复x值以避免插值问题
        unique_x, unique_indices = np.unique(x_sorted, return_index=True)
        x_sorted = x_sorted[unique_indices]
        y_sorted = y_sorted[unique_indices]

        if len(x_sorted) < 2:
            return 0.0, 0.0  # 不足点无法计算

        # 根据采样方法生成采样点
        if sampling_method == 'equidistant':
            # 等间距采样：在x方向上等间距采样
            x_interp = np.linspace(np.min(x_sorted), np.max(x_sorted), num_points)
            f = interp1d(x_sorted, y_sorted, kind='linear', fill_value='extrapolate')
            y_interp = f(x_interp)

        elif sampling_method == 'arc_length':
            # 等弧长采样：沿轮廓线弧长等间距采样
            dx = np.diff(x_sorted)
            dy = np.diff(y_sorted)
            ds = np.sqrt(dx ** 2 + dy ** 2)
            s = np.concatenate(([0], np.cumsum(ds)))  # 累积弧长
            total_arc_length = s[-1]

            if total_arc_length == 0:
                # 如果弧长为零，使用等间距作为后备
                x_interp = np.linspace(np.min(x_sorted), np.max(x_sorted), num_points)
                f = interp1d(x_sorted, y_sorted, kind='linear', fill_value='extrapolate')
                y_interp = f(x_interp)
            else:
                # 沿弧长等间距采样
                s_interp = np.linspace(0, total_arc_length, num_points)
                f_x = interp1d(s, x_sorted, kind='linear', fill_value='extrapolate')
                f_y = interp1d(s, y_sorted, kind='linear', fill_value='extrapolate')
                x_interp = f_x(s_interp)
                y_interp = f_y(s_interp)

        elif sampling_method == 'random':
            # 随机采样：从原始点中随机选择点
            if len(x_sorted) < num_points:
                indices = np.random.choice(len(x_sorted), size=num_points, replace=True)
            else:
                indices = np.random.choice(len(x_sorted), size=num_points, replace=False)
            x_interp = x_sorted[indices]
            y_interp = y_sorted[indices]
            # 排序以确保x递增
            sort_idx = np.argsort(x_interp)
            x_interp = x_interp[sort_idx]
            y_interp = y_interp[sort_idx]

        elif sampling_method == 'curvature_based':
            # 基于曲率的采样：在曲率高的区域增加点密度
            # 计算曲率（近似为二阶导数）
            dx = np.diff(x_sorted)
            dy = np.diff(y_sorted)
            if np.any(dx == 0):
                dx = np.where(dx == 0, 1e-5, dx)
            derivatives = dy / dx
            d2x = np.diff(x_sorted[:-1])
            d2y = np.diff(derivatives)
            if np.any(d2x == 0):
                d2x = np.where(d2x == 0, 1e-5, d2x)
            curvature = np.abs(d2y / d2x)

            # 归一化曲率并生成采样点权重
            curvature = np.append(curvature, 0)  # 对齐长度
            curvature = np.insert(curvature, 0, 0)  # 对齐长度
            weights = curvature / np.sum(curvature) if np.sum(curvature) > 0 else np.ones_like(x_sorted) / len(x_sorted)

            # 生成采样点：基于曲率权重
            indices = np.random.choice(len(x_sorted), size=num_points, p=weights, replace=True)
            x_interp = x_sorted[indices]
            y_interp = y_sorted[indices]
            # 排序采样点
            sort_idx = np.argsort(x_interp)
            x_interp = x_interp[sort_idx]
            y_interp = y_interp[sort_idx]

        else:
            raise ValueError(f"不支持的采样方法: {sampling_method}")

        # 计算Z2值
        dx = np.diff(x_interp)
        dy = np.diff(y_interp)
        if np.any(dx == 0):
            dx = np.where(dx == 0, 1e-5, dx)  # 避免除零
        derivatives_squared = (dy / dx) ** 2
        z2 = np.sqrt(np.mean(derivatives_squared))

        # 计算JRC
        jrc = 51.85 * (z2 ** 0.6) - 10.37
        return jrc, z2
    except Exception as e:
        print(f"JRC计算错误 ({sampling_method}): {str(e)}")
        return 0.0, 0.0


def create_visualization(enhanced_image: np.ndarray, fitted_curves: list, plot_output: str, image_width: int,
                         image_height: int) -> None:
    """创建对比可视化图像"""
    try:
        # 计算图像物理尺寸（毫米）
        image_width_mm = (image_width * MM_PER_INCH) / SINGLE_BOREHOLE_DPI
        image_height_mm = (image_height * MM_PER_INCH) / SINGLE_BOREHOLE_DPI

        fig, ax = plt.subplots(figsize=(16, 10))
        ax.imshow(enhanced_image, cmap='gray',
                  extent=[0, image_width_mm, image_height_mm, 0],
                  alpha=0.2, vmin=0, vmax=255)

        # 定义颜色序列
        fracture_colors = ['#FF6B6B', '#4ECDC4', '#45B7D1', '#FFD166', '#118AB2']

        # 绘制原始裂缝像素
        for i, (x_data, y_data, _) in enumerate(fitted_curves):
            ax.scatter(x_data, y_data, color=fracture_colors[i % len(fracture_colors)],
                       s=3, alpha=0.6, edgecolors='white', linewidth=0.3)

        # 绘制拟合曲线
        for _, (x_data, _, params) in enumerate(fitted_curves):
            x_plot = np.linspace(min(x_data), max(x_data), 200)
            y_plot = sine_func(x_plot, *params)
            ax.plot(x_plot, y_plot, color='black', linewidth=1.0)

        # 设置坐标轴
        ax.set_xlim(0, image_width_mm)
        ax.set_ylim(0, image_height_mm)
        ax.set_xlabel('周向距离 (mm)', fontsize=14, fontweight='bold')
        ax.set_ylabel('轴向深度 (mm)', fontsize=14, fontweight='bold')
        ax.set_title('裂隙识别与拟合结果对比', fontsize=16, fontweight='bold')

        # 优化刻度显示
        ax.grid(True, alpha=0.3, linestyle='--')
        ax.set_xticks(np.arange(0, image_width_mm + 1, 10))
        ax.set_yticks(np.arange(0, image_height_mm + 1, 50))

        # 添加图例
        from matplotlib.patches import Patch
        legend_elements = [
            Patch(facecolor='#FF6B6B', label='原始裂隙像素'),
            Patch(facecolor='black', label='拟合曲线')
        ]
        ax.legend(handles=legend_elements, bbox_to_anchor=(1.05, 1),
                  loc='upper left', fontsize=11, framealpha=0.9)

        # 保存图像
        plt.tight_layout()
        plt.savefig(plot_output, dpi=350, bbox_inches='tight', facecolor='white')
        plt.close()
    except Exception as e:
        print(f"可视化错误: {str(e)}")
        plt.close()


def process_image(image_path: str, output_csv: str = 'problem3-results.csv', plot_output: str = None) -> pd.DataFrame:
    """主处理函数：处理二值图像，进行裂隙聚类、正弦曲线拟合和JRC计算（支持多种采样方法）"""
    try:
        print(f"\n{'=' * 60}\n处理图像: {os.path.basename(image_path)}\n{'=' * 60}")

        # 读取图像
        image = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)
        if image is None:
            raise ValueError(f"无法读取图像 {image_path}")

        height, width = image.shape

        # 使用固定单孔DPI
        dpi = SINGLE_BOREHOLE_DPI
        print(f"使用单孔图像DPI: {dpi}")

        # 计算图像对应的物理尺寸
        image_width_mm = (width * MM_PER_INCH) / dpi
        image_height_mm = (height * MM_PER_INCH) / dpi
        print(f"图像物理尺寸: {image_width_mm:.2f}mm × {image_height_mm:.2f}mm")
        print(f"图像像素尺寸: {width} × {height}像素")
        print(f"钻孔周长: {BOREHOLE_CIRCUMFERENCE:.2f}mm")
        print(f"钻孔深度: {BOREHOLE_DEPTH}mm")

        # 图像增强
        enhanced_image = enhance_image(image)

        # 提取裂隙像素
        y_pixels, x_pixels = np.where(enhanced_image == 0)
        if len(x_pixels) == 0:
            print("警告：未检测到裂隙像素")
            return pd.DataFrame()

        print(f"检测到 {len(x_pixels)} 个裂隙像素")

        # 转换为物理单位（毫米） - 使用与文档1相同的转换方法
        x_coords_mm, y_coords_mm = pixel_to_mm_coordinates(
            x_pixels, y_pixels, width, height
        )

        # 连通区域分析
        mask = (enhanced_image == 0).astype(np.uint8)
        labeled_image = label(mask, connectivity=2)
        regions = regionprops(labeled_image)
        print(f"发现 {len(regions)} 个连通区域")

        # 处理每个区域
        fracture_list = []
        valid_count = 0

        for i, region in enumerate(regions):
            coords = region.coords
            y_pix, x_pix = coords[:, 0], coords[:, 1]

            # 转换为物理坐标 - 使用与文档1相同的转换方法
            x_phys, y_phys = pixel_to_mm_coordinates(
                x_pix, y_pix, width, height
            )
            region_coords = np.column_stack((x_phys, y_phys))

            is_valid, reason = validate_fracture(region_coords)
            if is_valid:
                fracture_list.append(region_coords)
                valid_count += 1
                print(f"区域 {i + 1}: ✓ 有效裂隙 ({reason})")
            else:
                print(f"区域 {i + 1}: ✗✗ 无效 ({reason})")

        print(f"\n有效裂隙数量: {valid_count}/{len(regions)}")

        # 正弦曲线拟合和JRC计算
        print("\n进行正弦曲线拟合和JRC计算...")
        results = []
        fitted_curves = []

        for idx, coords in enumerate(fracture_list):
            x_data = coords[:, 0]
            y_data = coords[:, 1]

            print(f"\n处理裂隙 {idx + 1}:")
            print(f"  像素数量: {len(x_data)}")
            print(f"  X范围: [{np.min(x_data):.2f}, {np.max(x_data):.2f}]mm")
            print(f"  Y范围: [{np.min(y_data):.2f}, {np.max(y_data):.2f}]mm")

            # 估计初始参数
            P_init = BOREHOLE_CIRCUMFERENCE
            C_init = np.median(y_data)
            R_init = (np.max(y_data) - np.min(y_data)) / 2
            beta_init = 0.0

            print(f"  初始参数: R={R_init:.2f}mm, P={P_init:.2f}mm, β={beta_init:.3f}rad, C={C_init:.2f}mm")

            # 使用修改后的拟合函数
            fit_result = robust_curve_fit(x_data, y_data, R_init, P_init, beta_init, C_init)
            R, P, beta, C, r_squared = fit_result

            print(f"  拟合参数: R={R:.2f}mm, P={P:.2f}mm, β={beta:.3f}rad, C={C:.2f}mm")
            print(f"  拟合质量: R² = {r_squared:.4f}")

            # 提取轮廓线并计算JRC（所有采样方法）
            region_mask = (labeled_image == regions[idx].label).astype(np.uint8)
            contour_points_pixel = extract_contour_skeleton(region_mask)
            if len(contour_points_pixel) == 0:
                print("  警告: 无法提取轮廓线，跳过JRC计算")
                jrc_results = {method: 0.0 for method in SAMPLING_METHODS}
                z2_results = {method: 0.0 for method in SAMPLING_METHODS}
            else:
                # 转换轮廓点坐标到物理单位 - 使用与文档1相同的转换方法
                x_contour, y_contour = pixel_to_mm_coordinates(
                    contour_points_pixel[:, 0], contour_points_pixel[:, 1], width, height
                )

                # 计算所有采样方法的JRC
                jrc_results = {}
                z2_results = {}
                for method in SAMPLING_METHODS:
                    jrc_value, z2_value = calculate_jrc(x_contour, y_contour, sampling_method=method,
                                                        num_points=DEFAULT_NUM_POINTS)
                    jrc_results[method] = jrc_value
                    z2_results[method] = z2_value
                    print(f"  {method}采样 JRC值: {jrc_value:.2f}, Z2值: {z2_value:.4f}")

            # 存储结果
            if r_squared >= 0.3:
                results.append({
                    '裂隙编号': idx + 1,
                    '振幅 R(mm)': round(R, 2),
                    '周期 P(mm)': round(P, 2),
                    '相位 β(rad)': round(beta, 4),
                    '中心线位置 C(mm)': round(C, 2),
                    'R_squared': round(r_squared, 4),
                    'JRC_等间距': round(jrc_results['equidistant'], 2),
                    'JRC_等弧长': round(jrc_results['arc_length'], 2),
                    'JRC_随机': round(jrc_results['random'], 2),
                    'JRC_基于曲率': round(jrc_results['curvature_based'], 2),
                    'Z2_等间距': round(z2_results['equidistant'], 4),
                    'Z2_等弧长': round(z2_results['arc_length'], 4),
                    'Z2_随机': round(z2_results['random'], 4),
                    'Z2_基于曲率': round(z2_results['curvature_based'], 4),
                    '像素数量': len(x_data)
                })
                fitted_curves.append((x_data, y_data, [R, P, beta, C]))
                print(f"  ✓ 拟合成功 (R²={r_squared:.4f})")
            else:
                print(f"  ✗✗ 拟合质量不足")

        print(f"\n拟合成功数量: {len(results)}/{len(fracture_list)}")

        # 输出结果
        if results:
            df = pd.DataFrame(results)
            output_columns = [
                '裂隙编号', '振幅 R(mm)', '周期 P(mm)', '相位 β(rad)', '中心线位置 C(mm)',
                'JRC_等间距', 'JRC_等弧长', 'JRC_随机', 'JRC_基于曲率',
                'Z2_等间距', 'Z2_等弧长', 'Z2_随机', 'Z2_基于曲率'
            ]
            df.to_csv(output_csv, index=False, columns=output_columns)
            print(f"结果已保存到 {output_csv}")

            # 显示结果摘要
            print(f"\n{'=' * 60}\n拟合结果摘要\n{'=' * 60}")
            for result in results:
                print(f"裂隙 {result['裂隙编号']}: R={result['振幅 R(mm)']}mm, P={result['周期 P(mm)']}mm, "
                      f"β={result['相位 β(rad)']}rad, C={result['中心线位置 C(mm)']}mm, "
                      f"JRC_等间距={result['JRC_等间距']}, JRC_等弧长={result['JRC_等弧长']}, "
                      f"JRC_随机={result['JRC_随机']}, JRC_基于曲率={result['JRC_基于曲率']}")
        else:
            print("无有效拟合结果")
            df = pd.DataFrame()

        # 可视化
        if plot_output and results:
            create_visualization(enhanced_image, fitted_curves, plot_output, width, height)
            print(f"可视化图已保存到 {plot_output}")

        return df

    except Exception as e:
        print(f"处理错误: {str(e)}")
        traceback.print_exc()
        return pd.DataFrame()


# 继续 process_all_images 函数
def process_all_images(input_dir: str, output_dir: str):
    """处理输入目录中的所有图像文件"""
    os.makedirs(output_dir, exist_ok=True)

    # 支持的图像扩展名
    image_extensions = ['.png', '.jpg', '.jpeg', '.bmp', '.tiff']
    image_files = [
        os.path.join(input_dir, f) for f in os.listdir(input_dir)
        if any(f.lower().endswith(ext) for ext in image_extensions)
    ]

    if not image_files:
        print(f"在目录 {input_dir} 中未找到图像文件")
        return

    print(f"找到 {len(image_files)} 个图像文件进行处理")

    # 处理每个图像
    all_results = []
    for image_path in image_files:
        base_name = os.path.splitext(os.path.basename(image_path))[0]
        output_csv = os.path.join(output_dir, f"{base_name}_results.csv")
        plot_output = os.path.join(output_dir, f"{base_name}_visualization.png")

        print(f"\n{'=' * 80}\n开始处理: {os.path.basename(image_path)}\n{'=' * 80}")
        results_df = process_image(image_path, output_csv, plot_output)

        if not results_df.empty:
            # 添加图像文件名到结果中
            results_df['图像文件名'] = os.path.basename(image_path)
            all_results.append(results_df)
            print(f"\n{'=' * 60}\n处理完成\n{'=' * 60}")
            print(results_df.to_string(index=False))
        else:
            print(f"\n处理完成，但未找到有效的正弦状裂隙")

        print(f"\n{'=' * 80}\n完成处理: {os.path.basename(image_path)}\n{'=' * 80}\n\n")

    # 保存所有结果的汇总
    if all_results:
        combined_df = pd.concat(all_results, ignore_index=True)
        summary_csv = os.path.join(output_dir, "all_results_summary.csv")
        combined_df.to_csv(summary_csv, index=False)
        print(f"所有结果汇总已保存到: {summary_csv}")

        # 输出采样方法比较统计
        print_sampling_method_comparison(combined_df)

    return all_results


def print_sampling_method_comparison(results_df: pd.DataFrame):
    """打印不同采样方法的比较统计"""
    print(f"\n{'=' * 60}\n采样方法比较分析\n{'=' * 60}")

    jrc_columns = ['JRC_等间距', 'JRC_等弧长', 'JRC_随机', 'JRC_基于曲率']

    for col in jrc_columns:
        if col in results_df.columns:
            valid_values = results_df[col][results_df[col] > 0]
            if len(valid_values) > 0:
                print(f"{col}:")
                print(f"  平均值: {valid_values.mean():.2f}")
                print(f"  标准差: {valid_values.std():.2f}")
                print(f"  最小值: {valid_values.min():.2f}")
                print(f"  最大值: {valid_values.max():.2f}")
                print(f"  有效样本数: {len(valid_values)}")
                print()


def create_sampling_comparison_plot(results_df: pd.DataFrame, output_path: str):
    """创建不同采样方法结果的比较图表"""
    try:
        jrc_columns = ['JRC_等间距', 'JRC_等弧长', 'JRC_随机', 'JRC_基于曲率']

        # 提取有效数据
        valid_data = []
        labels = []
        for col in jrc_columns:
            if col in results_df.columns:
                values = results_df[col][results_df[col] > 0].values
                if len(values) > 0:
                    valid_data.append(values)
                    labels.append(col.replace('JRC_', ''))

        if len(valid_data) < 2:
            return

        # 创建箱线图比较
        plt.figure(figsize=(12, 8))
        plt.boxplot(valid_data, labels=labels)
        plt.title('不同采样方法的JRC值分布比较', fontsize=16, fontweight='bold')
        plt.ylabel('JRC值', fontsize=14)
        plt.grid(True, alpha=0.3)
        plt.xticks(rotation=45)

        plt.tight_layout()
        plt.savefig(output_path, dpi=300, bbox_inches='tight')
        plt.close()
        print(f"采样方法比较图已保存到: {output_path}")

    except Exception as e:
        print(f"创建比较图错误: {str(e)}")


def analyze_sampling_method_differences(results_df: pd.DataFrame):
    """分析不同采样方法之间的差异"""
    print(f"\n{'=' * 60}\n采样方法差异分析\n{'=' * 60}")

    jrc_methods = {
        '等间距': 'JRC_等间距',
        '等弧长': 'JRC_等弧长',
        '随机': 'JRC_随机',
        '基于曲率': 'JRC_基于曲率'
    }

    # 计算两两方法之间的平均绝对差异
    methods = list(jrc_methods.keys())
    for i in range(len(methods)):
        for j in range(i + 1, len(methods)):
            col1 = jrc_methods[methods[i]]
            col2 = jrc_methods[methods[j]]

            if col1 in results_df.columns and col2 in results_df.columns:
                # 只比较两个方法都有有效值的情况
                mask = (results_df[col1] > 0) & (results_df[col2] > 0)
                if mask.sum() > 0:
                    diff = np.abs(results_df[col1][mask] - results_df[col2][mask])
                    mean_diff = diff.mean()
                    max_diff = diff.max()

                    print(f"{methods[i]} vs {methods[j]}:")
                    print(f"  平均绝对差异: {mean_diff:.2f}")
                    print(f"  最大绝对差异: {max_diff:.2f}")
                    print(f"  比较样本数: {mask.sum()}")
                    print()


# ==================== 扩展功能：裂隙粗糙度分析 ====================

def analyze_fracture_roughness(contour_points: np.ndarray) -> Dict[str, float]:
    """
    分析裂隙轮廓线的粗糙度特征
    参数:
        contour_points: 轮廓线坐标点数组
    返回:
        包含各种粗糙度指标的字典
    """
    if len(contour_points) < 2:
        return {}

    x = contour_points[:, 0]
    y = contour_points[:, 1]

    # 计算各种粗糙度指标
    results = {}

    # 1. 轮廓线长度与直线距离比
    dx = np.diff(x)
    dy = np.diff(y)
    contour_length = np.sum(np.sqrt(dx ** 2 + dy ** 2))
    straight_distance = np.sqrt((x[-1] - x[0]) ** 2 + (y[-1] - y[0]) ** 2)
    results['长度比'] = contour_length / straight_distance if straight_distance > 0 else 0

    # 2. 高度变化统计
    results['高度标准差'] = np.std(y)
    results['高度范围'] = np.max(y) - np.min(y)

    # 3. 斜率统计
    slopes = dy / np.where(dx == 0, 1e-5, dx)
    results['平均斜率'] = np.mean(np.abs(slopes))
    results['斜率标准差'] = np.std(slopes)

    return results


# ==================== 主程序 ====================
if __name__ == "__main__":
    print(f"{'=' * 60}\n钻孔裂隙图像处理系统启动\n{'=' * 60}")
    print(f"系统配置:")
    print(f"  输入目录: {os.path.abspath(INPUT_DIR)}")
    print(f"  输出目录: {os.path.abspath(OUTPUT_DIR)}")
    print(f"  单孔图像DPI: {SINGLE_BOREHOLE_DPI}")
    print(f"  钻孔直径: {BOREHOLE_DIAMETER}mm")
    print(f"  钻孔周长: {BOREHOLE_CIRCUMFERENCE:.2f}mm")
    print(f"  钻孔深度: {BOREHOLE_DEPTH}mm")
    print(f"  采样方法: {', '.join(SAMPLING_METHODS)}")
    print(f"  采样点数: {DEFAULT_NUM_POINTS}")
    print(f"{'=' * 60}")

    try:
        # 处理所有图像
        all_results = process_all_images(INPUT_DIR, OUTPUT_DIR)

        if all_results:
            # 创建采样方法比较图
            combined_df = pd.concat(all_results, ignore_index=True)
            comparison_plot = os.path.join(OUTPUT_DIR, "sampling_methods_comparison.png")
            create_sampling_comparison_plot(combined_df, comparison_plot)

            # 分析采样方法差异
            analyze_sampling_method_differences(combined_df)

            print(f"\n{'=' * 60}\n处理完成摘要\n{'=' * 60}")
            total_fractures = len(combined_df)
            print(f"总共处理的有效裂隙数量: {total_fractures}")

            # 统计各图像的裂隙数量
            if '图像文件名' in combined_df.columns:
                fractures_per_image = combined_df['图像文件名'].value_counts()
                print(f"\n各图像的裂隙数量:")
                for image, count in fractures_per_image.items():
                    print(f"  {image}: {count}个裂隙")

        print(f"\n{'=' * 60}\n所有处理完成!\n{'=' * 60}")
        print(f"结果保存在: {os.path.abspath(OUTPUT_DIR)}")

    except Exception as e:
        print(f"程序执行错误: {str(e)}")
        traceback.print_exc()
    finally:
        print(f"\n程序执行结束")